1. Import packages

library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6      ✔ purrr   0.3.5 
## ✔ tibble  3.1.8      ✔ dplyr   1.0.10
## ✔ tidyr   1.2.1      ✔ stringr 1.4.1 
## ✔ readr   2.1.3      ✔ forcats 0.5.2 
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
library(spData)
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
library(sf)
## Linking to GEOS 3.10.2, GDAL 3.4.2, PROJ 8.2.1; sf_use_s2() is TRUE
## New Packages
library(mapview) # new package that makes easy leaflet maps
library(foreach)
## 
## Attaching package: 'foreach'
## 
## The following objects are masked from 'package:purrr':
## 
##     accumulate, when
library(doParallel)
## Loading required package: iterators
## Loading required package: parallel
?registerDoParallel
registerDoParallel(4)
getDoParWorkers() # check registered cores
## [1] 4

2. Preparation for data

library(tidycensus)
racevars <- c(White = "P005003", 
              Black = "P005004", 
              Asian = "P005006", 
              Hispanic = "P004003")

options(tigris_use_cache = TRUE)
erie <- get_decennial(geography = "block", variables = racevars, 
                  state = "NY", county = "Erie County", geometry = TRUE,
                  summary_var = "P001001", cache_table=T) 
## Getting data from the 2010 decennial Census
## Using Census Summary File 1
buff=st_crop(erie,c(xmin=-78.9,xmax=-78.85,ymin=42.888,ymax=42.92))
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
buff$variable=as.factor((buff$variable))
variable_vec=as.array(unique(buff$variable))

rm(i)
## Warning in rm(i): object 'i' not found
points=foreach(i=1:4,.combine=rbind) %dopar%      
  {buff%>%filter(variable==unique(buff$variable)[i])%>%st_sample(size=.$value)%>%st_as_sf()%>%mutate(variable=unique(buff$variable)[i])}

3. Mpaview

mapview(points,zcol='variable',cex=1,alpha=0)